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Abstract 

The shear rate dependence of the shp length in thin polymer films confined between atomically 
flat surfaces is investigated by molecular dynamics simulations. The polymer melt is described 
by the bead-spring model of linear flexible chains. We found that at low shear rates the velocity 
proflles acquire a pronounced curvature near the wall and the absolute value of the negative slip 
length is approximately equal to thickness of the viscous interfacial layer. At higher shear rates, 
the velocity profiles become linear and the slip length increases rapidly as a function of shear rate. 
The gradual transition from no-slip to steady-state slip flow is associated with faster relaxation of 
the polymer chains near the wall evaluated from decay of the time autocorrelation function of the 
first normal mode. We also show that at high melt densities the friction coefficient at the interface 
between the polymer melt and the solid wall follows power law decay as a function of the slip 
velocity. At large slip velocities the friction coefficient is determined by the product of the surface 
induced peak in the structure factor, temperature and the contact density of the first fluid layer 
near the solid wall. 

PACS numbers: 68.08.-p, 83.80. Sg, 83.50.Rp, 47.61.-k, 83. ICRs 
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I. INTRODUCTION 



The rheology of complex fluids in thin films is important for theoretical and experimental 
studies of such common phenomena as friction, lubrication and wear Experimental 
measurements of the flow profiles and shear stresses on submicron scales might be subject 
to errors due to the possibility of liquid slip at the solid wall. An accurate prediction of flow, 
therefore, requires specification of a proper boundary condition. In the Navier model the 
interfacial shear rate and slip velocity are related via the proportionality coefficient, the so- 
called slip length. When the adjacent fluid layer slides past a solid wall with a finite velocity, 
the slip length is computed by linear extrapolation of the velocity profile near the interface 
to zero velocity [see Fig.[T](a)]. In the case when there is no relative velocity between fluid 
and solid at the interface, the formation of a lower viscosity boundary layer might result in 
the apparent slip length which is defined by the slope of the bulk velocity profile {2] [see 
Fig.[T](b)]. Experimental studies on pressure-driven flows in microchannels or thin film 
drainage using either the surface force apparatus (SFA) or atomic force microscope (AFM 
have demonstrated that the slip length depends on the nanoscale surface roughness 



surface wettability 



rate of shear js. 
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I2I . [13I . [1^ , and fluid structure 
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Several slip regimes exist due to progressive disentanglement of the anchoring chains in 



the shear flow of polymer melts 



13, 



18| . The difficulty in experimental determination of 



the velocity profiles and structure of complex fluids near interfaces leaves open important 
questions regarding the shear rate dependency of the slip length and the existence of a shear 
rate threshold for the boundary slip. 

In the last two decades, the boundary conditions at the interface between monatomic 



liquids and atomically 

n n r 

simulations [19!, 120, 121 



fat wa 
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Is were extensive 
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y studied by molecular dynamics (MD) 
j. The main factors affecting the slip are 



the energy of wall-fluid interaction and commensurability of the liquid and solid structures 
at the interface. At high wall-fluid energies, the first layer of fluid monomers becomes 
epitaxially locked to the solid substrate and the effective no-slip boundary plane is displaced 
into the fluid region. The absolute value of the negative slip length is approximately equal to 



the number of stacked monolayers between the effective and real boundary planes |19l . |21 |. 
At low surface energy, the first layer can slide with a finite velocity relative to the solid 
substrate under the shear stress from the bulk fluid [this situation is sometimes referred to as 
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a molecular or 'true' slip, e.g., see Fig.[T](a)]. The slip length is inversely proportional to the 
peak of the in-plane structure factor computed in the first fluid layer at the main reciprocal 



lattice vector 



21 



23, 



271 ]. The exact relation between the slip length and microscopic 



parameters of the liquid/solid interface, however, has not yet been established. 

The variation of the slip length with increasing shear rate in the flow of simple fluids pas t 
atomically smooth walls was first reported in the MD study by Thompson and Troian 29|. 
The nonlinear rate dependence of the slip length was well fitted by a power law function for 
different wall densities and weak wall-fluid interaction energies 29j. In a later study 27|], it 
was shown that the slip length is a linear function of shear rate at high wall-fluid interaction 
energies and, when the surface energy is reduced, the rate-dependence of the slip length can 
be well fitted by the power law function proposed in Ref. 29|]. It was also found that in a 
wide range of shear rates and wall-fluid interaction energies, the slip length is a function of a 
single variable which combines temperature of the first fluid layer, the contact density, and 
the peak value of the in-plane fluid structure factor evaluated at the main reciprocal lattice 
vector [27|. The results of previous MD studies of monatomic fluids lead to a conclusion 
that the boundary conditions for flows past smooth surfaces are either no-slip and rate- 
independent l2l| or described by the finite positive slip length which increases with shear 



rate 



28 



29l |. Except for flows of low density fluids [30[, there was no reported observation 
of a transition from no-slip to steady-state slip flow with increasing shear rate for simple 
fluids. 

The slip length in a flow of a polymer melt past a flat passive (non-adsorbing) surface is 
a ratio of the fluid viscosity to the friction coefficient which is determined by the interaction 
between fluid monomers and wall atoms 3l|]. In a recent MD study of unentangled melts 
confined between smooth surfaces 3^, it was found that the friction coefficient at the 



liquid/solid interface is nearly independent of the chain length beyond ten bead-spring units 
and, therefore, the molecular weight dependence of the slip length at low shear rates is mostly 



dominated by the melt viscosity. Depending on the strength o 



energy, the slip occurs either at the confining surfaces 32 



42 



33, 



the wa^ 
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431] or between the adsorbed layer and free polymer chains 



36 



1-fluid interaction 
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3511. A transition from 



stick to slip flow with increasinj 
hexadecane ^] and oligomers 



shear rate was reported in MD simulations of thin films of 
35| confined between strongly adsorbing surfaces. Despite 



extensive efforts in molecular simulations of thin polymer films, it still remains unclear what 
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system parameters (fluid density, pressure, chain length, surface energy) determine the onset 
of boundary slip at the interface between an unentangled melt and a flat surface. 

In a previous MD study [40'], the effect of shear rate on slip boundary conditions in thin 
polymer films confined between atomically smooth surfaces was investigated as a function 
of melt density. It was found that the slip length, extracted from the linear velocity profiles, 
passes through a local minimum at low shear rates and then increases rapidly at higher 
shear rates. This behavior was rationalized in terms of the friction coefficient (defined as a 
ratio of the shear-thinning viscosity to the slip length) , which undergoes a gradual transition 
from a nearly constant value to the power law decay as a function of the slip velocity. The 
functional form of the relation between the friction coefficient and the slip velocity in the 
crossover region can also be considered as a boundary condition for monatomic fluids 29|, |40 |. 
In addition, it was shown that in a wide range of shear rates and melt densities, the friction 
coefficient is determined by the product of the value of surface induced peak in the structure 
factor, temperature and the contact density of the first fluid layer near the solid wall 40] . 
Whether these conclusions remain valid at higher melt densities (when a viscous interfacial 
layer is formed near the confined surfaces) is one of the motivations of the present study. 

In this paper, we investigate the rate dependence of the slip length at the interface 
between an unentangled polymer melt and atomically flat walls using molecular dynamics 
simulations. We will show that in dense polymer films the velocity profiles are curved at low 
shear rates [shown schematically in Fig.[T](c) and Fig.[T](d)] due to highly viscous interfacial 
layer and the corresponding slip length is negative. When the viscosity of the interfacial 
layer is reduced at higher shear rates, the velocity profiles become linear and the slip length 
increases rapidly as a function of shear rate. The relaxation dynamics of polymer chains near 
the walls and in the bulk region is studied by analyzing the time autocorrelation function 
of the first normal mode at equilibrium and in shear flow. We will also show that the 
friction coefficient follows power law decay as a function of the slip velocity and its universal 
dependence on microscopic parameters of the liquid/solid interface holds only at large slip 
velocities. 

The rest of this paper is organized as follows. The details of molecular dynamics simula- 
tions and the equilibration procedure are described in the next section. The fluid velocity 
and density profiles, shear rate dependence of the slip length, and the analysis of the fluid 
structure and relaxation dynamics of polymer chains are presented in Section IIIII The 
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conclusions are given in the last section. 



II. MOLECULAR DYNAMICS SIMULATION MODEL 

I — I 

The simulation setup is similar to that described in the previous study [40] of a polymeric 
fluid undergoing planar shear flow between two atomically smooth walls. Figure[2] shows a 
snapshot of an unentangled polymer melt confined between solid walls. The total number of 
fluid monomers is kept the same as in the previous study {Nf = 6000) but the simulations 
are performed at higher fluid densities. The fluid monomers interact via the truncated 
Lennard- Jones (LJ) potential 

12 



VLj{r) = 4e 



r / Vr 



for r ^ rc = 2.5 cr, (1) 



where e and a are the energy and length scales of the fluid phase. The interaction between 
the wall atoms and fluid monomers is also modeled by the LJ potential with e^f = 0.9 e and 
o"wf = cr- The wall atoms are tethered about the sites of an fee lattice and do not interact 
with each other. 

The coarse-grained bead-spring model was used to represent an unentangled polymer melt 
with linear chains of = 20 monomers. In addition to the LJ potential any two neighboring 
monomers in the chain interact through the finitely extensible nonlinear elastic (FENE) 



potential 



44| 



V^ENE{r) = -^r^\n[l-ryr^], (2) 



with the standard parameters kg = 30 ea^"^ and Tq = 1.5 a [45]. A combination of the LJ and 
FENE potentials between neighboring monomers yields an effective spring potential, which 
is strong enough to prevent polymer chains from unphysical crossing each other or bond 
breaking even at the highest shear rates considered in the present study. 

The dynamics of fluid molecules and wall atoms was weakly coupled to a heat bath 



through a Langevin thermostat 46|]. The thermostat was applied only to the y component 
(perpendicular to the plane of shear) of the equations of motion for fluid monomers to avoid 
a bias in the shear flow direction [2l|]. The equations of motion for fluid monomers in the 
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X, y and z directions are given by 



myi + mVyi 



mxi 



mzi 




1 1 



(4) 



(5) 



(3) 



where the sum is taken over the neighboring fluid monomers and wall atoms within the cutoff 
radius rc = 2.5 o", F = 1.0 r^^ is the friction coefficient, and fi is a random uncorrelated force 
with zero mean and variance {fi{0)fj{t)) = 2mkBTT5{t)5ij determined from the ffuctuation- 
dissipation relation. The temperature of the Langevin thermostat is set to T = l.le/kB, 
where is the Boltzmann constant. The equations of motion were integrated using the 
fifth-order gear-predictor algorithm 4^] with a time step At = 0.002 r, where r= ^ma'^ je 
is the characteristic time of the LJ potential. The relatively small time step was chosen to 
resolve accurately the dynamics of fiuid molecules and wall atoms at the interface. 

The fiuid is confined by the fiat solid walls in the z direction (see Fig. [2]). Each wall 
consists of 576 LJ atoms arranged in two layers of an fee crystal with (111) plane parallel 
to the xy plane. The nearest-neighbor distance between fee lattice sites in the xy plane 
is (i = 1.0 0" and the wall density is pw = 1.40 cr~'^. The wall atoms were allowed to oscillate 
about their equilibrium lattice positions under the harmonic potential Vgp = | k with the 
spring stiffness k = 1200 e/cr^. The mean-square displacement of the wall atoms satisfies 
the Lindemann criterion for melting, i.e., {5u^)/d'^ < 0.023. In addition, the random force 
and the damping term were applied to all three components of the wall atom equations of 
motion, e.g., for the x component 



where the mass of the wall atoms is m^ = 10m, the friction coefficient is F = 1.0r~^ and 
the sum is taken over the fiuid monomers within the cutoff distance Tc = 2.5 cr. Periodic 
boundary conditions were applied in the x and y directions parallel to the walls. 

The polymer melt was initially equilibrated for about 5 x 10^ r at a constant normal 
pressure P = 0.5 ecr^^ applied on the upper wall while the lower wall was kept at rest. 
Then, the external pressure was gradually increased to a desired value (reported in Table[T]). 
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The distance between the walls was fixed after the system was additionally equilibrated 
at the constant pressure for about 2 x 10^ r. The MD simulations described below were 
performed at a constant density ensemble. The fiuid density, the corresponding pressure in 
the absence of shear flow, and the channel height are listed in TablelB 

The shear flow was generated by moving the upper wall with a constant velocity U in 
the X direction parallel to the stationary lower wall. Before the production run started, the 
flow was simulated for about 5 x lO^r for each value of the upper wall speed. The fluid 
velocity and density proflles were averaged within bins of thickness Az = 0.01 a for a time 
period up to 6 x lO^r at the lowest upper wall speed U = 0.001 cr/r. The small bin size was 
chosen to resolve accurately the velocity proflles and fluid structure near the walls. At the 
highest shear rates examined in this study, the Reynolds number did not exceed the value 
Re ^1.8 estimated from the maximum difference of the fluid velocities near the upper and 
lower walls, the shear viscosity, the fluid density and the channel height (see Table[T]). 



III. RESULTS 



A. Fluid density and velocity profiles 

Molecular simulations of polymer melts conflned by flat walls have demonstrated that 
the fluid structure consists of several discrete layers of monomers, the polymer chains are 
flattened near the walls, and the chain configurations remain bulklike in a region of several 



molecular diameters away from the walls [48|]. Qualitatively, these features can be observed 
in the snapshot of the polymer film shown in Fig.[2l Note that the first two layers of 
monomers near the walls are clearly distinguishable, the polymer chains in contact with 
the wall atoms are compressed towards the interface and their constitutive monomers are 
located within the discrete layers. 

The representative density profiles are shown in Fig.[3]for p = 1.0Aa^^ and p = 1.08cr^^ 
and two values of the upper wall speed. The density oscillations extend to a distance of 
about 5 — 6 cr from the walls and the profiles are uniform in the central bulk region. The 
magnitude of the first peak near the wall defines a contact density pc- The amplitude of the 
density oscillations is reduced with increasing upper wall speed. Notice that at low U in 
the case p = 1.08 cr~^ shown in Fig.[3](b), the local density minimum between the first and 
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second peaks is almost zero. It implies that the fluid monomers rarely jump between the 
first two layers, and it is, therefore, expected that the averaged velocity profiles will have 
relatively poor statistics in that region (see below). 

FigureH] shows the averaged velocity profiles in steady-state flow for the lowest fluid den- 
sity (p = 1.04(j~^) considered in this study. The velocity profiles remain linear throughout 
the channel except for a larger slope inside the first fluid layer. A finite slip velocity is 
noticeable even at very low upper wall speed U= 0.01 a/r. The data are noisy because 
the averaged velocity component in the x direction is much smaller than the thermal fluid 
velocity = kBT/m. The velocity profile for U= 0.01 a/r bends slightly near the walls, 
which implies that the interfacial viscosity is higher than the fluid bulk viscosity. A small 
curvature in the bulk region of the velocity profile at U= Aa/r might be related to the 



40[ |. The normalized slip velocity 



nonuniform heating up of the fluid at high shear rates 
increases with increasing upper wall speed. 

The averaged velocity profiles at the higher melt density (p = 1.08cr~^) are reported in 
Fig.[5l At small values of the upper wall speed, the slip velocity of the first fluid layer is 
barely noticeable and the profiles are highly curved near the walls due to the presence of 
the viscous interfacial layer with thickness of about 2 — 3a. The statistical fluctuations are 
relatively large near the walls because of the pronounced density layering [e.g., see Fig.[3](b)]. 
The effective no-slip boundary plane is located inside the fluid domain at a distance of about 
2 a from the wall. With increasing upper wall speed, the fluid velocity profiles become linear, 
the no-shp boundary plane is displaced out of the fluid region, and the slip velocity increases. 

The normalized velocity profiles for the highest fluid density p = 1.11 a^'^ are plotted in 
Fig.O At the lowest upper wall speed U = 0.01 a/r, the first four monolayers stick to the 
walls but the velocity profiles remain linear in the middle of the channel. Due to the low 
probability of finding chain segments in between the monolayers the statistical uncertainties 
are much larger near the walls than in the bulk region. We found that the shape of the flow 
profiles depends on how the system was prepared. In one case, the upper wall speed was 
increased from zero to U = 0.01 cr/r following the equilibration procedure described in the 
previous section. The averaged velocity profile is marked by the lower black curve in Fig. [61 
In the other case, the upper wall speed was first increased toU = 0.5 cr/r, and then, after the 
equilibration period of about 10^ r, gradually reduced to U = 0.01 a/r. The corresponding 
velocity profile is shown by the upper red curve in Fig.O Note that the width of the flowing 
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regions is nearly the same in both cases but the thickness of the immobile interfacial layers 
varies from about four to five molecular diameters. We also found that at higher upper 
wall speeds (0.1 < f/r/cr < 0.25) the first fluid layer slides with a finite velocity past the 
substrate and the weak oscillations in the velocity profiles near the walls correlate with the 
fluid density layering (see inset in Fig.E]). Finally, as shown in the inset of Fig.El the velocity 
profiles become linear throughout the channel at higher upper wall speeds {U > 0.3 a/r) 
and the slip velocity increases monotonically with increasing U. 

The slope of the linear part of the velocity profiles in the bulk region of the channel (6 a 
away from the walls) was used to compute the shear rate, shear stress and slip length. We 
define the slip length as a location of the point where linearly extrapolated velocity profile 
vanishes. Negative value of the slip length implies that the effective no-slip boundary plane 
is displaced into the bulk fluid domain. The velocity of the first fluid layer with respect to 
the lower stationary wall was computed as follows 

V, = rV,{z)p{z)dz I r p{z)dz, (7) 

where the limits of integration {zq = —14.50 a and Zi = — 13.85 cr) were determined from 
the width of the first peak in the density profile. Note that even if the flow profiles are 
linear throughout the channel, the velocity of the first layer Vi is slightly larger than the slip 
velocity computed from the Navier relation Vg = jLg (by about 0.5 7 units of cr/r). This 
is because the slip length is defined with respect to the reference plane 0.5 a away from the 
inner fee lattice plane and the first fluid layer is located approximately a away from the fee 
plane (see Fig. [3]). 



B. Shear rate dependence of the melt viscosity and slip length 

We first estimate the polymer melt viscosity which is defined as a ratio of shear stress to 



shear rate, i.e., axz = pj- The shear stress was computed using the Kirkwood relation 49 1 
in the bulk region {6 a away from the confining walls), where the fluid structure is uniform 
and the velocity profiles are linear even at low shear rates. The viscosity is plotted in 
Fig-El as a function of shear rate for the indicated values of the polymer density. The 
Newtonian regime is observed only in a narrow range of shear rates, and it is followed by 
the crossover to the shear-thinning behavior, which occurs at higher shear rates when the 
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melt density is reduced. The dashed hne in Fig. [7] corresponds to the power law decay with 
the exponent —0.37 reported for the lower density polymer melts (0.86 ^ pcr^ ^ 1.02) at 
high shear rates 40]. The data presented in Fig. [7] indicate that the melt viscosity in the 
shear-thinning regime decreases faster at higher fluid densities. The statistical errors due to 
thermal fluctuations are relatively large at low shear rates. 

The variation of the slip length as a function of shear rate is presented in Fig. [8] for all 
melt densities considered. As expected from the shape of the velocity profiles described in 
the previous section, the slip length at low shear rates is negative (except for p = 1.04a~^) 
and its magnitude is approximately equal to the thickness of the viscous interfacial layer. 
With increasing shear rate, the velocity profiles become linear, implying that the local 
viscosity of the boundary layer is reduced, and the slip length increases rapidly. At the 
highest melt density p = l.llcr~^ and low shear rates, the thickness of the interfacial layer 
depends on the equilibration procedure and the slip length cannot be uniquely defined. The 
uncertainty in the slip length related to the thickness of the immobile interfacial layer and 
the slope of the velocity profile in the bulk region is about 2 a. At higher upper wall speeds 
(0.1 "^U T / a < 0.25) and p= 1.11 a"'^, the first fluid layer is sliding with a finite velocity and 
the slip length is relatively large. Multivalued slip lengths were also reported in shear flow 
of simple fluids past smooth surfaces with high wall-fluid interaction energy 2l|. Finally, we 
comment that the data presented in Fig. [8] cannot be well fitted by the power law function 
proposed in Ref. [291]. 

In our simulations, the upper wall speed was varied so that the slip velocity remains 
less than about the fluid thermal velocity. We performed test runs at higher upper wall 
speeds {U > 4cr/r) and observed a different regime, where the the slip length becomes a 
nonmonotonic function of shear rate. It was also recently shown that the slip length at the 
interface between short chain polymers and smooth thermal walls approaches a constant 
value at high shear rates 50|, ISjJ . In the present study, the behavior of the slip length at 
very large slip velocities and high shear rates was not examined in detail. 

The rate-dependent boundary conditions can be reformulated in terms of the wall shear 
stress and slip velocity. In the steady-state shear flow, the stress across any plane parallel to 
the confining walls is the same, and, therefore, the shear stress computed in the bulk region 
is equal to the wall shear stress. The slip velocity of the first fluid layer was calculated by 
averaging the velocity profile over the width of the first density peak using Eq. ([71). Note 
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that the shp velocity computed from the Navier relation {Vs = "yLs) is smaller than Vi or 
can be even negative if the velocity profiles are curved near the interface [e.g., see Fig.[T](c) 
and Fig.[T](d)]. The shear stress averaged in the bulk region is plotted in Fig.lHlas a function 
of the slip velocity of the first fiuid layer for the indicated polymer melt densities. The 
shear stress increases rapidly at small slip velocities and grows steadily at large Vi. At the 
highest melt density p=l.lla~^, the time-averaged shear stress is discontinuous at small 
slip velocities. 

The nonlinear relation between the shear stress and slip velocity shown in Fig.[9]can be 
used to determine the friction coefficient per unit area axz = kVi as a function of Vi . We 
comment that the ratio of the shear viscosity to the slip length cannot be used to compute 
the friction coefficient at the liquid/solid interface when the velocity profiles are curved 
near the .urface^nd the dip length i. extracted from the bulk part of the prohle. In the 
previous study [40], the simulations were performed at lower polymer melt densities and 
the velocity profiles remained linear at all shear rates examined. In the range of densities 
(0.86 ^ pcr^ ^ 1-02), the friction coefficient {k = fJ./Ls) as a function of the slip velocity was 
well described by the following equation 

k/k* = [i + ivjv:)']-'-'', (8) 

where k* and V* are the normalization parameters [4^. In the present study, the friction 
coefficient is plotted as a function of the slip velocity in Fig.[TOl The data can be well fitted 
by the empirical formula Eq. (jH]) at lower melt densities p ^ 1.06 a^^ (not shown). At higher 
melt densities p> 1.08 0"^"^, the plateau regime, where the friction coefficient is independent 
of the slip velocity, is absent and the slope of the power law decay is slightly smaller than 
—0.7 (see the dashed line in Fig. fTU]) . 



C. Friction coefficient and fluid structure in the first layer 



The connection between friction at the liquid/solid interface and fiuid structure induced 



by the periodic sur 



polymer melts 



32 



ace potential was established for monatomic fiuids 



21 



40| , adsorbed monolayers 5J, |55| and adsorbed polymer layers 



271, 



53|, 



The 



11 



in-plane structure factor in the first fiuid layer near the solid substrate is defined as 

where rj = {xj,yj) is the position vector of the j-th monomer and Ni is the number of 
monomers within the layer j21^]. Typically, the structure factor exhibits a circular ridge at 
the wavevector |k| ^ 271/ a due to short range ordering of the fiuid monomers. In addition, 
at sufficiently high wall-fiuid interaction energy, several sharp peaks appear in the structure 
factor at the reciprocal lattice vectors of the crystal wall. It is well established that the mag- 
nitude of the peak at the first reciprocal lattice vector in the shear fiow direction correlates 
well with the friction coefficient at the interface between a simple LJ liquid and a solid wall 



composed out of periodically arranged LJ atoms [2l|, |23|, l27| . 



In the present study, the first reciprocal lattice vector in the shear fiow direction 
Gi = (7.23 (7^^, 0) is slightly displaced from the wavevector |k| ^ 27r/cr. The averaged struc- 
ture factor ^(Gi) is plotted in Fig.[TT](a) as a function of the slip velocity of the first fiuid 
layer for the polymer density p = 1.08 cr~'^. The error bars are relatively large at small slip 
velocities due to the slow relaxation dynamics of the polymer chains in the interfacial layer 
(see next section). At higher slip velocities, fiuid monomers spend less time in the min- 
ima of the periodic surface potential, and, as shown in Fig.[TT](a), the magnitude of the 
surface-induced peak in the structure factor decreases logarithmically with increasing Vi. 
Furthermore, as expected from the density profiles shown in Fig. [3] for different upper wall 
speeds, the contact density is reduced at higher slip velocities [see Fig.[TT](b)]. Similarly to 
the definition of the slip velocity given by Eq. ([7]), the temperature of the first fiuid layer 
was computed as follows 

Ti = /" T{z)p{z)dz I I p{z)dz, (10) 

J ZO 'J ZO 

where T{z) is the local kinetic temperature and the limits of integration (2:0 = —14.50 cr and 
Zi = — 13.85 0") were determined from the width of the first density peak. The variation of 
the monolayer temperature as a function of the slip velocity is presented in Fig.[TT](c). At 
small slip velocities Vi < 0.03 cr/r, the temperature is equal to value T = 1.1 e/ks set by the 
Langevin thermostat. With increasing upper wall speed, the temperature of the first fiuid 
layer gradually rises up to T ^1.6e/kB at the highest slip velocity Vi ~ 1.3 cr/r reported 
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in Fig.[TT](c). At shear rates 7 > 0.01 r~^, the temperature profiles across the channel T{z) 
become nonuniform and the heating up is larger near the interfaces j4o|. 

The dependence of the friction coefficient on the structure factor, contact density and 



temperature of the first fiuid layer was studied in the recent paper [401 at lower polymer 
melt densities (0.86 ^ pcr^ ^ 1.02). Except for the densities p = 1.00cr~^ and 1.02a^^ at 
low shear rates, the data for the friction coefficient {k = fi/Ls) were found to collapse on 
master curves when plotted as a function of either Ti/ [S (Gi) pc] or S{0)/ [S{Gi) pd- In 



both cases, the data cou 
and —1.15, respectively 



d be well fitted by a power law function with the exponents —0.9 
40l ]. At higher polymer densities (1.04 ^ pa^ ^ 1.11) considered 
in the present study, the inverse friction coefficient is plotted in Figures [T^] and dSl as a 
function of the combined variables Ti/ [S (Gi) pc] and S{0)/ [S{Gi) pc], respectively. The 
collapse of the data holds at small values of the friction coefficient k < Aera'"^ and the 
surface-induced peak in the structure factor ^(Gi) < 14. Thus, the simulations at higher 
melt densities provide an upper bound for the friction coefficient {k ^ 4eTa~^), below which 
the data are described by a single master curve [40|. These results support the conclusion 
from the previous studies that at the interface between an atomically smooth solid wall and 
a simple |27] or polymeric [40] fluid, the friction coefficient is determined by a combination 
of parameters evaluated in the flrst fluid layer. 



D. Relaxation dynamics of polymer chains 

Equilibrium molecular dynamics studies of polymer melts conflned between attractive 
walls have shown that the relaxation of polymer chains slows down near the interfaces and 
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becomes bulk-like at distances of about two radiuses of gyration away from the walls 
In order to probe the relaxation dynamics of polymer chains in shear flow we evaluated the 
autocorrelation function of the normal modes in the direction perpendicular to the plane 
of shear. The y component of the normal coordinates for a discrete polymer chain of 
monomers is given by 

YpW = ^E^^w-«^^^r^' (11) 

i=l 

where i/i is the component of the position vector of the i-th monomer in the chain and 



p = 0, 1, A^ — 1 is the mode number 



59( 1 . The longest relaxation time of a polymer chain 
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is associated with the first mode p = 1. The normahzed time autocorrelation function for 
the first normal mode is computed as follows 

Ci(t) = (Yi(t)-Yi(0))/(Yi(0)-Yi(0)). (12) 

The relaxation dynamics in heterogeneous systems is usually described by the stretched 
exponential (or Kohlrausch- Williams- Watts) function Ci{t) = exp[— (t/ri)^]. The time in- 
tegral of the stretched exponential defines the characteristic decay time r* = tiT{P~^)//3, 
where F is the gamma function and (3 is the stretched exponential coefficient. 

The autocorrelation function, Eq. (fT2|) . was computed separately in the bulk region {6 a 
away from the solid walls) and inside the interfacial layers (within 3a from the walls). 
During the simulation, the normal coordinate and the position of the center of mass of each 
polymer chain were calculated every 100 r, and the autocorrelation function was updated 
only if the position of the center of mass was inside either bulk or interfacial regions. The 
correlation function was averaged for at least 6 x lO^r at low shear rates to resolve very slow 
dynamics of the polymer chains near the interfaces. 

Figureim shows the time autocorrelation function computed at equilibrium conditions 
(i.e., U = 0) for the indicated polymer densities. As expected, the relaxation of the polymer 
chains in the bulk and near the walls is slower at higher melt densities. The inverse bulk 
relaxation time, estimated roughly from Ci{t) = 1/e in Fig.[T3](a), correlates well with the 
onset of shear thinning of the fiuid viscosity reported in Fig. [71 For each density, the decay 
of the correlation function is slower for the chains in the interfacial layer than in the bulk of 
the channel. The difference in relaxation times is especially evident at higher melt densities, 
p ^ 1.09 cr~^, where only the early stage of relaxation is reported in Fig.[14](b). At the 
lowest polymer density p= 1.04 cr~^, the decay time of the interfacial chains is about 2 — 3 
times larger than in the bulk, which agrees with the conclusion drawn from the shape of 
the velocity profile for U = 0.01 a/r in Fig.lHthat the interfacial viscosity is higher than the 
bulk value. 

The infiuence of the shear fiow on the time autocorrelation function is demonstrated in 
Fig.fTSl for the polymer density p = 1.08cr~^. At low upper wall speeds U < 0.1 a/r, the 
relaxation dynamics near the interface is much slower than in the bulk region, implying the 
formation of a highly viscous interfacial layer. The spatial variation of the shear viscosity 
correlates with the nonlinearity in the velocity profiles shown in Fig.O When U= 1.0 a/r, 
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the decay time of the correlation function is nearly the same across the channel (see the 
dashed curves in Fig. fTSl) and the velocity profile is linear (shown for U = 0.5 a/r in Fig.[5]). 
At higher upper wall speeds, the polymer chains relax slightly faster near the walls than in 
the bulk, possibly because of the heating up of the fluid near the walls. 

Finally, the relaxation time of the polymer chains near the walls and the slip length 
are summarized in Fig. [16] as a function of shear stress and polymer density. The leftmost 
points of the curves shown in Fig. [16] correspond to the upper wall speed U= 0.01 cr/r. At 
low shear stress, the relaxation time varies widely from r* ^ 6.2 x 10^ r at p= 1.04 o"^^ to 
r* ^ 3.3 X 10^ r at p= 1.11 a~^, indicating the presence of a highly viscous boundary layer 
at higher melt densities. With increasing shear stress, the decay time decreases, and, when 
U > 1.5 cr/r, the relaxation of the polymer chains near the walls becomes even slightly faster 
than in the bulk region (see also Fig. [T5|) . The dependence the inverse relaxation time on 
the shear stress exhibits qualitatively similar behavior to the slip length, which supports the 
conclusion (drawn from the shape of the velocity profiles) that the transition to slip flow is 
associated with shear- melting of the boundary layer. The inset in Fig. [16] shows the variation 
of the stretched exponential coefficient as a function shear stress. The data are scattered 
at low shear stress due to the slow relaxation of the polymer chains, and the decay of the 
autocorrelation function becomes nearly exponential at higher shear stress. 



IV. CONCLUSIONS 



In this paper, the rate-dependence of the slip length at the interface between a dense 
polymer melt and weakly attractive smooth walls was studied using molecular dynamics 
simulations. The melt was modeled as a collection of linear chain polymers (A^ = 20). It 
was shown that at low shear rates the velocity profiles are curved near the wall due to 
the formation of a highly viscous interfacial layer and the effective slip length is negative 
and almost rate-independent. With increasing upper wall speed, the gradual transition to 
steady-state slip fiow is associated with shear-melting of the interfacial layer. The relaxation 
dynamics of polymer chains in shear fiow was analyzed by evaluating the decay of time 
autocorrelation function of the first normal mode in the vorticity direction. We found that 
the rate behavior of the slip length correlates well with the inverse relaxation time of the 
polymer chains in the interfacial layer. 
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The rate-dependent slip boundary conditions were also reformulated in terms of the 
friction coefficient at the polymer/wall interface and slip velocity of the first fiuid layer. In 
agreement with the results of the previous study [40], we found that the friction coefficient 
at lower melt densities undergoes a transition from a constant value to the power law decay 
as a function of the slip velocity. At higher melt densities the friction coefficient decays as 
a power law function in a wide range of slip velocities. When the magnitude of the surface 
induced peak in the fluid structure factor is below a certain value, the friction coefficient 
is determined by a combination of parameters (structure factor, temperature, and contact 
density) of the flrst fluid layer near the solid wall. 
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Figures 
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FIG. 1: (Color online) A schematic diagram of a steady-state shear flow over a flat solid wall, (a) 
The slip velocity and the slope of the linear velocity profile are related via Vs = jLg, where Lg is 
the slip length, (b) Apparent slip is associated with the lower viscosity boundary layer, (c) The 
velocity profile with a downward curvature due to higher viscosity boundary layer is described by 
the negative slip length, (d) A combination of the 'true' slip at the liquid/solid interface and the 
curvature of the velocity profile. 
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FIG. 2: (Color online) A snapshot of fluid monomers (open blue circles) and wall atoms (filled 
gray circles) positions projected on the xz plane. The bottom wall is at rest and the top wall is 
moving with a constant velocity U in the x direction. Each fluid monomer belongs to a polymer 
chain. Five polymer chains are marked by solid lines. The fluid monomer density is p = 1.08(T~^ 
and the top wall speed is U = O.lcr/r. 
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FIG. 3: (Color online) Averaged monomer density profiles near the stationary lower wall with 
e^i/e = 0.9. The uniform fluid density away from the walls \s p= 1.04 o""^ (a) and p= 1.08 cr^^ (b). 
The upper wall velocities U are tabulated in the inset. The left vertical axis indicates the location 
of the fee lattice plane at z = — 15.24 o". The dashed line at z= — 14.74(7 denotes the reference 
plane for computing the slip length. 
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FIG. 4: (Color online) Average normalized velocity profiles for the indicated values of the upper 
wall speed and the fluid density p = 1.04 o"^^. The vertical axes coincide with the position of the 
fee lattice planes at z= — 15.24 a and z = 8.77a. The dashed lines denote liquid/solid interfaces 
at z = — 14.74: a and z = 8.27 a. The inset shows an enlarged view of the velocity profiles near the 
stationary lower wall. 
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FIG. 5: (Color online) Averaged velocity profiles for the indicated upper wall speeds U (in units 
of (t/t) and the fluid density p= 1.08 o"^^. The vertical axes indicate the location of the fee lattice 
planes at 2;= — 15.24(T and 2: = 7.86(7. The dashed lines denote the position of the liquid/solid 
interfaces 0.5 cr away from the fee planes. The region near the lower wall is enlarged in the inset. 



24 




-12 



-8 



-4 

z/o 



p = 


1.1 la"^ 


= u 


= 0.01 a/x 


— u 


= 0.125 a/x 


— u 


= 0.3 a/x 




= 1.0 a/x 


1 1 1 1 





4 



FIG. 6: (Color online) Average normalized velocity profiles for the upper wall speed U = 0.01 cr/r 
(see text for description). The vertical dashed lines indicate liquid/solid interfaces at 2;= — 14.74(j 
and z = 6.81 a. The inset shows averaged velocity profiles for the upper wall speeds U = 0.125 fi/r 
(dashed violet curve), U = 0.3 u/r (continuous blue curve), and U = 1.0 cr/r (dash-dotted green 
line) . 
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FIG. 7: (Color online) Shear rate dependence of the bulk viscosity /x/ero"^'^ for the indicated 
values of the polymer density. The dashed line with a slope —0.37 is shown for reference. Solid 
curves are a guide for the eye. 
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FIG. 8: (Color online) Slip length Lgja as a function of shear rate for the indicated values of the 
polymer melt density. At low shear rates and p = \.\\a^^ the slip length is multivalued. Solid 
curves are drawn to guide the eye. 
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FIG. 9: (Color online) Shear stress (in units of ea~^) averaged in the bulk region as a function 
of the slip velocity of the first fluid layer for the indicated values of the polymer melt density. Solid 
curves are a guide for the eye. 
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FIG. 10: (Color online) Log-log plot of the friction coefficient k = (Jxz/Vi (in units of ero""^) as 
a function of the slip velocity [computed from Eq. ([7])] for the tabulated values of the polymer 
density. The dashed line with a slope —0.7 is plotted for reference. The same data as in Fig.[9l 
Solid curves are drawn as a guide for the eye. 
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FIG. 11: (Color online) Structure factor 5(Gi) evaluated at the first reciprocal lattice vector 
Gi = (7.23(7^^,0) in the shear flow direction (a), contact density near the stationary lower wall 
(b), and temperature (c) of the first fluid layer as a function of the slip velocity Vi (in units of 
cr/r). The solid blue line y = 6.93 — 7.73 ln(x) represents the best fit to the data. Solid curves are 
a guide for the eye. 
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FIG. 12: (Color online) Log-log plot of the inverse friction coefficient k^^ = Vijoxz (in units of 
a'^/er) as a function of Tiks/s [S{Gi) pccr^]~^ evaluated in the first fluid layer for the indicated 
polymer melt densities. The solid line with a slope 0.9 is plotted as a reference. 
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FIG. 13: (Color online) Log-log plot of the inverse friction coefficient k^^ = Vijoxz (in units of 
cr^/er) as a function of the variable >S'(0)/ [S'(Gi) pc o'^\ computed in the first fluid layer. The fluid 
densities are listed in the inset. The dashed line with a slope 1.15 is shown for reference. 
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FIG. 14: (Color online) Normalized time autocorrelation function of the first normal mode at 
equilibrium for polymer chains in the bulk region (a) and near the walls (b) for the indicated 
values of the fluid density. 
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FIG. 15: (Color online) Time autocorrelation function of the first normal mode for polymer chains 
in the bulk region (a) and near the walls (b) and the fluid density />= l.OSo"^'^. The upper wall 
speed from right to left is C/ = 0, 0.025, 0.1, 0.2, 0.5, 1.0, 2.0, 3.0 and 4.0 (in units of <t/t). 
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FIG. 16: (Color online) Inverse relaxation time of the polymer chains near the walls (top) and the 
slip length (bottom) as a function of shear stress axz (in units of ea~^) for the indicated values of 
the melt density. The inset shows the stretched exponential coefficient for the same data. 
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Tables 
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TABLE I: The fluid density p and pressure P at equilibrium (U = 0). The density is defined as a 
ratio of the total number of fiuid monomers {Nf = 6000) to the volume 20.86 a x 12.04 a x (h — a), 
where h is the distance between the fee lattice planes in contact with the fiuid. 



p (units of cr~^) 


1.04 


1.06 


1.08 


1.09 


1.11 


P (units of £a~^) 


6.0 


7.0 


8.0 


9.0 


10.0 


h (units of a) 


24.01 


23.50 


23.10 


22.83 


22.54 
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